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Strain engineering has emerged as a powerful tool to modify the optical and electronic properties of 
two-dimensional crystals. Here we perform a systematic study of strained semiconducting transition 
metal dichalcogenides. The effect of strain is considered within a full Slater-Koster tight-binding 
model, which provides us with the band structure in the whole Brillouin zone. From this, we derive 
an effective low-energy model valid around the K point of the BZ, which includes terms up to second 
order in momentum and strain. For a generic profile of strain, we show that the solutions for this 
model can be expressed in terms of the harmonic oscillator and double quantum well models, for 
the valence and conduction bands respectively. We further study the shift of the position of the 
electron and hole band edges due to uniform strain. Finally, we discuss the importance of spin-strain 
coupling in these 2D semiconducting materials. 

PACS numbers: 73.63.-b,68.65.-k,71.70.Fk,71.70.Di,11.15.Wx 


I. INTRODUCTION 

The outstanding stretchability of the new families of 
2D crystals makes them excellent candidates for their 
use in strain engineering [1]. This opens the possibility 
to fabricate nanodevices in which the optical and elec¬ 
tronic properties are tunable by controlled application 
of external strain. Single layer of semiconductor transi¬ 
tion metal dichalcogenides (TMDs) with the form MX 2 
(where M=Mo, W is a transition metal and X=S, she 
is a chalcogen atom) can sustain large amounts of strain 
before rupture of the membrane. For these materials, 
a direct-to-indirect band gap transition is expected for 
uniaxial/biaxial tensile strain of the order of ~ 2 — 3% 
[2, 3], and a semiconducting-to-metal transition has been 
predicted for 10 — 15% of tensile biaxial strain [4, 5]. 

Like graphene, low-energy excitations of insulating 
MX 2 are mainly localized close to the two inequivalent 
points K and K', also denoted as “valleys”, paving the 
way for the possibility of valleytronics, namely convey 
the information in the valley degree of freedom. The pe¬ 
culiarities of these materials suggest also possible ways to 
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manipulate the valley-bit. The effect of strain on stan¬ 
dard silicon semiconductor physics is known to lead to an 
enhancement of the electron and hole mobilities, and to 
the valley-degeneracy breaking [6]. However, silicon de¬ 
vices cannot sustain strains larger than ~ 1.5%, whereas 
single layer TMDs support strength deformations higher 
than 10% [7, 8]. The strong spin-orbit coupling (SOC) 
indeed yields a different spin-polarization of the valence 
band. Therefore, several degrees of freedom are strongly 
entangled in TMDs [9, 10]. Tuning the spin-orbit cou¬ 
pling of mechanical deformation has been explored in 
conventional GaAs based semiconductors and quantum 
wells where a linear strain dependence is found in this 
coupling [11-13]. Indeed very recently, a coupling be¬ 
tween single electron spins and the motion of mechanical 
resonators based on crystal strain has been reported ex¬ 
perimentally [14]. Therefore, controlling and tailoring 
their properties, at the applicant as well as at the theo¬ 
retical level, represents thus the current challenge for a 
wide community of scientists. 

In the last years, graphene has become the natural 
platform to test strain engineering physics. One of the 
main theoretical study of deformed graphene-based ma¬ 
terial was done by Kane and Mele [15] who used a tight- 
binding model to study the effect of long wavelength 
deformations on the low-energy electronic structure of 



2 


carbon nanotubes. They showed that the effects of the 
tubule geometrical features and symmetry on its elec¬ 
tronic structure are included through an effective vector 
potential. Such gauge held has been also predicted by 
Suzuura and Ando [16] in the context of electron-phonon 
scattering in carbon nanotubes, and a group symmetry 
based survey has been done by Manes [17]. For specihc 
prohles of strain, it was predicted theoretically [18] and 
then proved experimentally [19-21] that pseudo Landau 
level quantization corresponding to strong effective mag¬ 
netic helds can be realized in graphene. Moreover, this 
kind of pseudomagnetic field is also observed in an arti- 
hcial molecular graphene assembled by atomic manipu¬ 
lation of carbon monoxide molecules over a conventional 
two-dimensional electron system on a copper surface [22] . 

Strain engineering methods have been applied to other 
2D crystals, and recently the possibility to tune the band 
gap with strain has been experimentally proven for M 0 S 2 
[23-26] and WS 2 [27-29]. Moreover, spatially modu¬ 
lated biaxial tensile strain has been applied to single 
layer M 0 S 2 , leading to the realization of an optoelectronic 
crystal consisting of artificial atoms, due to the spatial 
modulation of the band gap in the sample [30]. Piezo¬ 
electricity and piezoresistivity effects have been recently 
reported for single layer and multi-layer M 0 S 2 [31, 32]. 
Therefore, there is a need for a deeper understanding of 
the effect of external non-uniform strain on the physical 
properties of semiconductor TMDs. 

Here, we theoretically investigate the effect of strain 
on the electronic structure of a monolayer MA 2 . Our 
main focus will be to study the effect of inhomogeneous 
strain on the low-energy physics of the system. We 
start by considering a Slater-Koster tight-binding model 
which contains the relevant orbital character in the va¬ 
lence and conduction band, originated from dxy 

and dx^-y 2 orbitals of the M metal atom, and Px, Py and 
Pz orbitals of the chalcogen atom X [10, 33]. Strain is 
considered in this model by means of a modification of 
the corresponding hopping terms [25] . From that model, 
we use the Lowdin partitioning method [34] to obtain 
an analytical two-bands k • p Hamiltonian valid in the 
vicinity of the K points of the BZ. This model differs on 
the previously used Dirac-like Hamiltonian [35, 36] in the 
fact that we include terms up to second order in momen¬ 
tum and strain, which are needed to capture some of the 
main features of the MX 2 electronic band dispersion at 
low energies. One of the consequences is that strain leads 
to the appearance of not only one pseudo gauge field in 
the theory, but also to the existence of several pseudo 
vector potentials that couples in the relevant terms in 
the low-energy theory for MX 2 , and which are absent in 
the well known case of strained graphene [37, 38]. 

We further apply our model Hamiltonian to the case of 
arc-shaped deformation, and find that the corresponding 
solutions are well described in terms of the harmonic os¬ 
cillator and double quantum well models, for the valence 
and conduction bands, respectively. Finally, we study 
the strain-induced change of the valence and conduction 


band edges, as well as the coupling between spins and 
strain on this family of TMDs. 

The paper is organized as follows. In Sec. H we con¬ 
sider the effect of strain within a full Slater-Koster tight- 
binding model. From it, a low-energy k • p Hamiltonian 
is derived in Sec. HI, and the impact of pseudo-vector 
fields on the electronic properties is discussed in Sec. IV. 
In Sec. V we do an analytical analysis of the low-energy 
electronic spectrum of single layer TMDs using single 
band pictures. In Secs. VI and VH we discuss the effect 
of strain on the position of the valence and conduction 
band edges, and spin-orbit coupling, respectively. Our 
main results are summarized in Sec. VHI. 


II. TIGHT-BINDING MODEL FOR STRAINED 
TMDS 

Monolayer M 0 S 2 is a direct band gap semiconductor 
with the gap placed at the K and K’ points of the hexag¬ 
onal BZ. Ab-initio calculations show that there are two 
additional secondary extrema: a local maximum of the 
valence band at the F point, and a local minimum of 
the conduction band at approximately at the Q point, 
midway between F and K point [39]. These features, 
which are not relevant to the main optical properties of 
the system, might play an important role in transport 
properties [40, 41]. The low-energy physics of monolayer 
M 0 S 2 around the K and K’ points was first described 
by a simple massive Dirac Hamiltonian [35]. More accu¬ 
rate approximations have been developed later, as tight- 
binding methods [33, 42-44] and k • p approximations 
[42, 45] which goes beyond the massive Dirac model, and 
account for the presence of trigonal warping and diag¬ 
onal quadratic terms in momentum. In this section, we 
describe the TB theory that will be used as starting point 
to consider strain effects on the electronic band structure 
of M 0 S 2 . 

The main features of the band structure of monolayer 
M 0 S 2 in the whole BZ are well captured by the Slater- 
Koster TB model of Ref. [33] which includes 11 bands 
corresponding to the d orbitals of the metal atom and to 
the p orbitals of the chalcogen atoms. Remarkably, the 
relevant physics of monolayer M 0 S 2 around the gap is 
covered by a smaller subspace, which can be obtained by 
performing an appropriate unitary transformation that 
transform the p orbitals of the top and bottom S lay¬ 
ers into their symmetric and antisymmetric combinations 
with respect to the z-axis. For the single-layer case, the 
resulting 11-band model can be decoupled in 6 bands 
with even symmetry with respect to z —> —z inversion, 
and 5 bands with odd symmetry. Low-energy excitations 
belong exclusively to the first block, so that we will dis¬ 
regard the other states with different symmetry. Local 
spin-orbit interaction can be as well included in a suit¬ 
able way [10]. Diagonal terms oc LzSz appear here to 
be dominant, so that in good approximation each spin 
sector can be dealt with separately [10]. Using the com- 


3 



FIG. 1. (Color online) A top view schematic of monolayer 
M 0 S 2 lattice structure. Blue (Orange) circles indicate Mo 
(S) atoms. The nearest neighbors {Si) and the next nearest 
neighbors (oi) vector have been shown in the figure. 


Crystal Fields 

eo 

-1.094 


£2 

-1.512 


ep 

-3.560 


tz 

-6.886 

Intralayer Mo-S 

^pdcr 

3.689 


^pdTT 

-1.241 

Intralayer Mo-Mo 

^ddcr 

-0.895 


VddTT 

0.252 


Vdds 

0.228 

Intralayer S-S 

^ppa 

1.225 


VppTT 

-0.467 


TABLE I. Slater-Koster tight-binding parameters for single¬ 
layer M 0 S 2 . All terms are in units of eV. 

provide with the different contributions in Appendix A. 
An appropriate set of Slater-Koster parameters for M 0 S 2 
is given in Table I. 


pact notation of Ref. [33], we can consider the reduced 

Hilbert space: A. Hamiltonian in strained lattice 


(^ 32 ^ —da,2_y2, dxy:Px: Py ^ Pz') (^) 

where the S and A superscripts of the p-orbitals re¬ 
fer to symmetric and antisymmetric combinations pf = 
lj'/^{pl +Pi) and pf = lly/2{p\ — p\), the index i runs 
over the spatial directions i = x,y, z, and the labels t and 
b refer to the top and bottom sulfur planes, respectively. 
A top view of the crystal lattice of M 0 S 2 is sketched in 
Fig. 1. The TB Hamiltonian defined by the base (1), 
including the local spin-orbit coupling can be expressed 
in the real space as 

H = 'y ] -f y ] + H.C.], (2) 


where c| ^ creates an electron in the unit cell i in the 
atomic orbital labelled hy p = 1,... ,6 belonging to the 
Hilbert space defined in Eq. (1). The Hamiltonian ac¬ 
quires a more compact form once written in the fc-space: 

TT _ ( Hmm Hmx\ 

\Hmx^ Hxx J ’ 

F^MM = eM + 2 ^ t“^cos(k-aj), 

Hxx = ^x + ‘i ^ tf^cos(k-ai), 

2^1,2,3 

Hmx = 1 : (3) 

1=1,2,3 

where the nearest {8i) and the next nearest (a^) neigh¬ 
bour vectors are shown in Fig. 1. All the hopping 
terms have been evaluated within a Slater-Koster 

scheme [10, 25, 33, 46]. For the sake of simplicity, we 


The use of a Slater-Koster tight-binding approach is 
particularly convenient when lattice deformations, like 
strain, are considered. Neglecting as a first approxima¬ 
tion the corrections to the local atomic potentials due 
to lattice deformation [47, 48] the effect of strain is here 
driven by the dependence of the tight-binding parameters 
of the two-center energy integral elements which depend 
on the interatomic distance. The effect of strain is thus 
considered in the simplest way by varying the interatomic 
bond lengths as a result of the applied strain. At the lin¬ 
ear order, the modified hopping terms in the presence of 
strain can be written as 


j 1 Aij^ 




( 4 ) 


where jrF] is the distance in the equilibrium positions 
between two atoms labelled by {i,p) and {j,iz), jr^-j 
the distance in the presence of strain, and = 

—dlntij_^i^(r)/dln(r)j,r=|r°.| is the dimensionless bond- 
resolved local electron-phonon coupling. For practical 
purposes, we have [r?-] = ^JljVl a for the M-X bond, 
and jrFj = a for the in-plane M-M and X-X bonds. 

A microscopic evaluation of the electron-lattice cou¬ 
pling parameters 

is in principle possible by means of an accurate analy¬ 
sis based on the direct comparison between DFT and 
tight-binding calculations. Along this line, for instance, 
the electron-phonon coupling associated with the differ¬ 
ent interlayer hopping in multilayer graphene were es¬ 
timated in Ref. [49]. This task turns however to be 
formidable in transition metal dichalcogenides because 
of the large number of orbitals/bands and because of the 
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lack of a Fermi surface that can be used as a reference. In 
the absence of any theoretical and experimental estima¬ 
tion for the electron-phonon coupling, we use the Wills- 
Harrison argument [50] namely oc 

where is the absolute value of the angular momen¬ 
tum of the orbital /i, and li, is the absolute value of the 
angular momentum of the orbital v. Following this ap¬ 
proach we assume that Aij^x-x = 3, Aij^x-M = 4, and 
Aij^M-M = 5, for the X-X pp, for X-M pd, and for 
the M-M dd hybridizations, respectively. The applica¬ 
tion of strain transforms the vector Tq, which separates 
two lattice sites connected with electron hopping, into 
r « ro -I- rg • Vu. Note that, in the above transforma¬ 
tion, we are considering only the acoustic part of the 
displacement vector, which has been shown to be a good 
approximation in the long wavelength region of interest 
here [16]. 

In general, we can write Vu = e + u>, where e and 
uj are the strain and rotation tensors, respectively [48]. 
The strain tensor of a two-dimensional system is given 
by the symmetric tensor 



with components that include both in-plane and out-of- 
plane displacements 

1 f dui duj \ 1 duz 

= 2 8 ^ 3^7 <**> 

where r = {x,y) and u = {ux,Uy,Uz) are the position 
and displacement vectors, respectively. The rotation ten¬ 
sor UJ accounts for local rotations in the system. It is 
an antisymmetric tensor defined by 2ujxy = —2ujyx = 
{duy/dx — dux/dy). For a homogeneous strain the rota¬ 
tion tensor is zero and we can assume thatr = ro-(l-l-£). 
On the other hand, for an inhomogeneous strain with the 
local rotation we must use r = rg • (1 -I- Vu) [48]. Ex¬ 
plicit expressions for the atomic separation as modified 
by non-uniform strain are given in Appendix B. 

III. LOW-ENERGY MODEL OF STRAINED 
TMDS 

Hamiltonian (2), which includes explicitly the hy¬ 
bridization between the metal and the chalcogen atoms, 
represents the appropriate starting point for a compelling 
derivation of an effective low-energy model in the pres¬ 
ence of strain. For this purpose, we perform a Taylor ex¬ 
pansion in momentum and in the strain fields, followed 
by a canonical projection onto the two (conduction and 
valence) low-energy bands. From the technical point of 
view, in order to obtain an effective 2x2 (4x4 including 
spin) model Hamiltonian, we use the Lowdin partitioning 
method [34]. Details about the derivation are provided in 
Appendix B. Similar to the carbon nanotube and to the 
graphene cases [16, 51-53] we first set the momentum co¬ 
ordinates on the relative valley (K-point of the Brillouin 


Ao = —O.lleF Ao = 69 meV 
A = 1.82eF A = —81 meV 
Ao = — 17meF A' = —2 meV 



eV 


eV 

meV 


meV 

0 + 

15.99 

Qi 

15.92 

«r 

-61 


-5.7 


-3.07 

0.2 

-1.36 

a“+ 

3.2 

«2 

0.02 


-0.17 

Qg 

0.0 

a^+ 

3.4 


0.01 


hi 

V2 

V3 

Vi 

75 

0.002 

-56.551 

1.635 

1.362 

8.180 


TABLE II. Microscopical parameters of the spinful two- 
band low-energy model. The upper table describes strain- 
independent Hamiltonian parameters where to = 2.34 eV, 
a = —0.01 and j3 = —1.54; the middle table the Hamilto¬ 
nian parameters related to the strain through a scalar poten¬ 
tial [Eq. (8)]; the lower table the Hamiltonian parameters rji 
related to the strain-induced coupling to the pseudo-vector 
potentials A^. 

zone), and we derive hence a strain-dependent Hamilto¬ 
nian which includes the effect of hopping integrals modifi¬ 
cation caused by the deformation. The strain-dependent 
Hamiltonian around K-point, up to second order in strain 
and momentum, can be written as H = Hq -|- Hgo, where 

Ho = + V + tgOg (^q -b • (Tr 

+ (|q + -b [q -b IrAgl^/Scr^) , 

Hso = I -b (5A -b Og (^[q -b ^tA41^Ao 

-b|q-b ^rAgpA'cr^^ I TS. (7) 

Here e is the elementary charge, rug is the free elec¬ 
tron mass, (Jx = {xUx,(Jy) are Pauli matrices in the 
2x2 “band” space, and s = ± and r = ± are spin 
and valley indexes, respectively. Finally og = a/y/S and 
q = is the relative momentum with respect to 

the K point. The parameters Ag, A, Ag, A, Ag, A', tg, a 
and /3 are strain-independent and they can be obtained 
directly from the Slater-Koster parameters of the origi¬ 
nal Hamiltonian (2), and they are given in Table H for 
the case of M 0 S 2 . A detailed derivation of the numeri¬ 
cal values of all the parameters of the low-energy model 
in terms of the original tight-binding parameters can be 
found in Appendix B. It is useful to notice that the mass 
asymmetry parameter, a, and topological term, /3, are re¬ 
lated to general physical properties of the band structure, 
like effective mass and energy gap, through the relations 
a = mg/m+ and (3 = mo/m- — dmgu^/(A — A_), where 
V = tgOg/S, m± = mcmv/(rnv ± me), 2A± = Ag ± A. 
In addition, me and are the effective masses of the 
conduction and valence band, and A+ and A_ are the 
spin-orbit coupling of the conduction and valence bands. 
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respectively [42]. 

The presence of a finite strain induces in the Hamilto¬ 
nian (7) many different terms. The most straightforward 
are the diagonal ones, i.e. a scalar potential, which con¬ 
tains a spin independent part, V — diag[I?+,D_], and a 
spin-dependent contribution, SX = diag[i5A+, dA_]. The 
explicit expressions of D± and <5A± read: 

D± = af [Alp + ujly) + afV^, 
dX±=af\A\^ + afiV + culy) + afV^. ( 8 ) 

Note that the strain fields appear in Eqs. (7) and (8) 
only through the representative variables A = Sxx—Syy — 
^‘X‘Cxyi ^ ^xx 4“ ^yyi ^cid ujxy (^dtiyIdx duxjdy^l^j. 
The numerical values of all are also reported in Ta¬ 
ble. II. It should be noticed that the quantitative use of 
the second order terms in the scalar potential {V and 
(5A) should be done with some cares. Because, here, we 
considered the linear approximation to include deforma¬ 
tion in the bond lengths. However, these terms would be 
negligible for small deformation. 

In addition to the above discussed diagonal terms, it is 
interesting to underline the appearance in (7) of five dif¬ 
ferent fictitious gauge fields defined as A, where 

Ax = (/i/eao)Re[Al] and Ay = {h/eao)lm\A\. The cou¬ 
pling constants rji are evaluated from the values of the 
initial Slater-Koster parameters, and their specific value 
for the case of single-layer M 0 S 2 are reported in Table H. 
Note that, due to the small value of rji, the off-diagonal 
pseudo vector potential (Ai) results to be very weak as 
compared to the diagonal ones. The opposite happens for 
the well known cases of mono- and bilayer graphene, for 
which the off-diagonal terms are the dominant compo¬ 
nents of the strain dependent Hamiltonian [16, 37]. The 
weakness of Ai in M 0 S 2 might stem from the large energy 
gap as compared with graphene, which is a semi-metal 
with no gap. 


IV. STRAINED TMDS AS A 
MULTI-PSEUDO-VECTOR FIELD SYSTEM 

The dependence of the electronic/transport/optical 
properties of TMDs triggers the biggest interest towards 
realistic applications for strain-engineering in these ma¬ 
terials and hence many theoretical and experimental se¬ 
tups have been proposed. Most of them employ the de¬ 
pendence of the magnitude of the optical or transport 
gap. On the conceptual basis, such proposals are thus 
related to the strain modulation of the scalar potentials 
D±, (5A±. Interesting enough, there is, on the other hand, 
off-diagonal terms that can be described in terms of the 
pseudo-vector potential. The concept of pseudo gauge 
fields, for instance, has been widely discussed in the con¬ 
text of strained graphene [16, 37, 51, 52] and it provides 
the possibility to induce extremely large effective pseudo- 
magnetic fields [19]. Such pseudomagnetic fields are thus 
reflected in the onset of flat bands (Landau levels) in the 


energy spectrum of deformed systems, as observed ex¬ 
perimentally in strained graphene samples [19, 20, 22]. 
A similar framework has been discussed in TMDs by 
Cazalilla et aZ. [36]. However, the Hamiltonian in Eq. (7), 
appropriate for realistic modelling of monolayer TMDs, 
shows profound differences in regards to this simplistic 
picture since at least three pseudo vector potentials (Ai, 
A 2 , and A 3 ) are induced by strain, even in the sim¬ 
plest spinless case. It should be stressed that, due to 
this multi-pseudo-vector field structure, these fields can¬ 
not be referred as gauge fields, since the contemporary 
presence of three fields cannot be described by a simple 
phase shift in the wave-function after doing a transfor¬ 
mation like A — > A -I- VA. In other words, the effect of 
these pseudo vector fields can not be eliminated by the 
counter-acting presence of a real magnetic field. 

The complexity of these multi pseudovector field struc¬ 
tures gives rise to qualitatively new physics, which is not 
present in graphene-like systems and in the previous anal¬ 
yses for TMDs based on only one vector potential. The 
rich phenomenology of this structure will be thus the ob¬ 
ject of investigation of the present Section. 

An interesting feature in graphene under strain is the 
possibility of tailoring the pseudo vector potential in 
snch a way to mimic the effects of an effective mag¬ 
netic field. Under these conditions, (pseudo-) Landau 
levels are expected to appear, reflected in flat bands 
in the electronic band structure. This scenario has 
been theoretically predicted [18] and experimentally ver¬ 
ified in graphene [19] and similar predictions have been 
prompted out in TMDs [36]. Things are actnally more 
complex in a realistic modelization of the TMDs, due to 
the presence of many pseudo-vector potentials as we will 
discuss later. 

In order to focus on the possible occurrence of pseudo 
Landau levels (PLLs), we neglect in the following the role 
of the scalar potentials, and we consider only the leading 
term of the spin-orbit coupling in the absence of strain. 
Only three pseudo vector potentials will appear, Aj with 
j = 1, 2, 3. The first standard step to address this issue is 
to introduce the total canonical momentum fields iTj — 
{hq eAj) = {hq erjjA), where tvj = nf -\- zttJ. The 
fundamental thing to be underlined here is that the fields 
TTj are not orthogonal, but fulfill the following relations: 

[tTj, TTj] = -ieh{rij - rii){dx -k idy){Ax + iAy) 

ki, 7rj] = -ieh{r]j - ? 7 i)V • A - eh{rij -k x A)^. 

( 9 ) 

It is interesting to notice that when all vector poten¬ 
tials have the same couplings rji = t], then there is just 
one gauge field, which leads to the standard algebra for 
the canonical momentum associated with a real magnetic 
field. Therefore, this kind of solution corresponds to the 
real Landau levels of the system in the presence of a true 
magnetic field. 

However, the commutation relations in (9) for the more 
general case of strained TMDs, imply that such oper- 
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ators do not commute and obey a more complex alge¬ 
bra. Finding an analytical solution for the PLLs results 
thus a formidable task due to the non-orthogonality of 
the theory [54] which is one of the main consequences 
of the presence of multi-pseudo-vector fields. Neverthe¬ 
less, it is instructive to consider the symmetric gauge, 
A = ^{—y,x), which is the one associated with the ex¬ 
perimentally relevant case of trigonal deformation of the 
lattice [18, 19]. In this case one can show that 

[di,oJ]=5'ij , [oi,%]=0 , [a|',a]]=0 (10) 


where Stj = {rji + r]j)/2^/\ri, 


tVj\ 


and we have intro¬ 


duced the creation operators dj = 


y/Whin. 


TT,, where 


Ib = is the magnetic length in terms of the 

pseudo-magnetic field B = |V x A|. The overlap matrix 
of this case, Sij ^ implies that the bosonic opera¬ 
tors Qi are in general non-orthogonal, and [hi,hj] yf 0 
for i ^ j. In order to give a solution for this problem, 
one needs to redefine Fock space in a non-orthogonal ba¬ 
sis [54]. In Sec. V we will provide with a perturbative 
analytical solution using a single band model. 


allow for an analytical determination of the energy 
bands and of the pseudo Landau levels. We have 
thus solved the problem numerically. In order to re¬ 
veal the relevant physics associated with the pseudo¬ 
vector potentials and with possible pseudo-Landau lev¬ 
els, we choose an inhomogeneous strain profile corre¬ 
sponding to a constant pseudo magnetic field for each 
vector potential. An arc-shape deformation (sketched 
in Fig. 2), which corresponds to a displacement profile 
{ux,Uy) = {xy/R,—x^/2R) (i? being here the arc radius), 
is known to be one of the simplest and efficient candi¬ 
dates [55, 56]. Within this context the resulting gauge 
field in the Landau gauge reads A = -^{y/R,0), corre¬ 
sponding to a three different constant pseudo magnetic 
field Bi = rjih/eaQR, associated to the pseudo vector po¬ 
tentials Ai = r]iA. Neglecting the weak contribution of 
the rotation tensor, Hamiltonian (7) can be written in 
first-quantization as 


A. Energy spectrum in the Landau gauge: 
arc-shape deformation 

Besides this specific case, however, the complex multi¬ 
vector potential structure of the Hamiltonian does not 


+ toaoiQx + - toaod. 


H = 


aoR' 


^toaoiqx + Vi^) + toaody V-{y) - ^(a - f3)d^ 


where 


Tx / \ (^0 Aqs) dz (A + As) 7 -, / \ / \ 

V±{y) = —--Z + D±{y) +SX±{y)s+ - - a [ Qx + m 

2 Amo 


( 11 ) 




The solution of the above Hamiltonian leads to a set of 
coupled differential equations that we solve numerically 
for hard-wall boundary conditions. Details can be found 
in Appendix D. 

The resulting dispersion relations under different strain 
conditions are shown in Fig. 3 and (4) where, in order 
to analyze the possible presence of pseudo-Landau levels, 
the scalar potentials {D± and 6X±) have been neglected. 
In comparison, we show the band dispersion of the un¬ 
strained case in Fig. (3)a. Apart from the parabolicity 
of the valence and conduction bands, one can observe a 
crossing of the edge modes. This is expected due to the 
nonzero Chern number associated to each flavour (spin 
or valley) in our model (i.e. 2C = sign(A) — sign(/3)) [57]. 
In other words, as long as the valley index is a good 


quantum number (for instance in zigzag termination and 
hard-well boundary cases), and the valley-Chern number 
is integer valued, then metallic edge modes are expected 
to exist in the gap. Although this problem is beyond 
the scope of the present paper, the nature of these edge 
states will be discussed in Sec. V C using the numerical 
results of TB model. Here we mention that the existence 
of metallic edge modes in this kind of systems depend on 
the edge potential and atomic termination. [58, 59] 

Fig. 3b shows the band dispersions in the previously 
discussed toy model where the coupling of all the three 
strain induced fields is the same {r]i = 77 ), which corre¬ 
spond to the case of a real magnetic field applied to the 
sample. Remarkably, the first Landau level in the va¬ 
lence band evolves in an electron-like edge mode. This 
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FIG. 2. (Color online) Top view of an arc-shaped MX 2 with 
R = 4Ly in which blue and red lattice points indicates M and 
X atoms, respectively. 
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FIG. 3. (Color online) (a) Energy dispersion of unstrained 
nanoribbon case with Ly = 50ao. (b) Energy dispersion a 
strained nanoribbon for Ly = 50ao, rit = rj = 1 and R = 
0.5Ly. Noticeably, the first Landau level in the valence band 
evolves in an electron-like edge mode. 


feature is consistent with the results of the full six-band 
tight-binding model [46]. 

The presence of flat bands (pseudo-Landau levels) can 
be speculated from these results. However, as we are go¬ 
ing to show, the actual relevance of pseudo-Landau levels 
appear doubtful in the more realistic case of arc-shaped 



qxao 

EIG. 4. (Color online) Energy dispersion of arc-shaped case 
with scalar potential and rotation tensor neglected. The dis¬ 
persion relation show equidistant parabolic valence subbands, 
and a set of conduction subbands that present band cross¬ 
ing for the low-energy states and parabolic dispersion above 
some specific energy. The size of the energy gap and level 
spacing between the subbands is strongly dependent on the 
strain strength. We set (&)Ly = 50ao,R — 2.5Ly and (b) 
Ly = 50ao, R = 1.8Ly. 


tension. In Fig. 4 we set the parameters rji with the nu¬ 
merical values listed in Table II, which correspond to the 
realistic full tight-binding dispersion. Panels (a) and (b) 
show the band dispersion of a given spin flavour for two 
different magnitudes of strain, parametrized in terms of 
two values of the arc radius i?. The resulting dispersion 
show almost equidistant parabolic valence subbands, and 
a set of conduction subbands that present band cross¬ 
ing for the low-energy states and parabolic dispersion 
above some specific energy. It should be noticed also 
that the size of the energy gap and level spacing between 
the subbands (in both valence and conduction sectors) is 
strongly dependent on the strain strength. In the follow¬ 
ing section we will show how these peculiar features, for 
both valence and conduction bands, can be understood 
in terms of a harmonic oscillator and inverted harmonic 
oscillator physics, respectively. 





















In summary, comparing these results in strained M 0 S 2 
(as a representative case of single layer TMD), with 
those of strained graphene [37] reveals three main dif¬ 
ferences between these two systems: i) No obvious sim¬ 
ilarity between pseudo vector potentials and real mag¬ 
netic field, as it happens in strained graphene, ii) Strong 
particle-hole asymmetry in the strained M 0 S 2 energy 
dispersion, as compared to the symmetric spectrum in 
strained graphene, iii) Absence of flat PLL in M 0 S 2 
for Landau’s gauge (i.e. arc-shape) deformation in con¬ 
trast with the appearance of flat bands in the arc-shaped 
graphene [55, 56]. 


V. SINGLE-BAND (EFFECTIVE MASS) 
MODEL FOR STRAINED TMDS 


The two-band Hamiltonian (7) can be perturbatively 
decomposed into two individual one-band Hamiltonians 
for the conduction and the valence bands at small q and 
strain. Those single-band Hamiltonian which is valid 
around the K point of the BZ, have the following analyti¬ 
cal expressions (see Appendix C for a detailed derivation) 


H±=E±^ D± 

± 7 |q+|Air} 


r I ® A I 

--{a|q-f -A 2 I 

477iQ fi 


±^|q+ 


h 


(13) 


iTl I Ao±A (tpap) 

2 2 ^ 'n (A-x 


7 = 


where E± = —^ ^ '/i (A-.A-)i| ’ 

4moU^/(A — sA_), s is the spin index and -|-(—) indi¬ 
cates conduction (valence) band, respectively. Having 
neglected a and /3 terms, we could get the model Hamil¬ 
tonian used in Ref. [36] to study a pseudomagnetic field in 
a monolayer M 0 S 2 , which reveals (in the absence of SOC) 
symmetric PLL in the conduction and valence bands. 
Notice that the single band model (13), considers the ef¬ 
fective mass asymmetry (a) and momentum dependent 
mass term (/3), leading to different model solutions for 
the two cases. 

The Hamiltonian from (13) can be easily deduced for 
the Landau’s gauge (arc-shape) deformation and by ne¬ 
glecting scalar potential contribution D± we have 


pseudo vector potentials, depending on the sign of {a ± 
P ± ^)wf' in the single-band models. One is a harmonic 
oscillator (HO) physics where (Q!±/3±7)wf > 0 and the 
other is a double quantum well (DQW) physics where 
{a ± fi ± < 0. In the case of M 0 S 2 which is 

addressing here, after plugging the numerical value of 
the model parameters in (14), we obtain {wf,wf) ~ 
(-32.2,-24.0) X (aoR)-"^, (u>^,w^) ~ (0.06,0.12) x 
(aoR), (ic^, wff) ~ (4.05, —3.57), a + fi + ^ = 3.92 and 
a — fi — j = —3.94 for the up component of spin index. 
Therefore, the model for M 0 S 2 leads to a DQW and HO 
for the conduction and valence bands, respectively. 

We would like to emphasize that the single band 
Hamiltonian (14) provides with a general model that can 
be applied to other families of strained semiconductor 
2D crystals. For the sake of completeness, we describe 
the generalities of the model. First of all, we point out 
that a + fi + ^ > 0 and a — /I — 7 < 0, since these signs 
originate from the positive and negative masses of elec¬ 
trons and holes, respectively. These relations imply that 
0 < I a] < fi + j which imposes some restriction over 
the two-band model parameters just based on the sign of 
effective masses. 

The sign of however, might change depending on 
the microscopic properties of the considered system. Fur¬ 
thermore, if the effective masses are the same for both 
bands (i.e. a = 0) then wf = —wf. In this case, there 
are two possibilities, namely > 0 {wi < 0) which 
leads to HO (DQW) solutions for both the valence and 
conduction bands. Therefore, an asymmetry in the ef¬ 
fective masses is necessary to have two different physics 
(understood as HO or DQW energy spectra) in the elec¬ 
tron and hole bands. 

Finally, we briefly discuss the case of trigonal defor¬ 
mation {ux,Uy) = ^ (^xy, ^ ^, where uq have the 

units of inverse length and quantifies the strength of the 
applied strain. This strain profile has been widely dis¬ 
cussed in the context of graphene [18, 19] and recently 
it has been considered in TMDs [36]. Such a deforma¬ 
tion is properly described by the symmetric gauge, i.e. 
A = ^^(y, —x), which leads to the single band Hamil¬ 
tonian 


H±=E± 


4m( 


■ [(a ± /3 ± 7 )g^ +w^(y- w^q^f + w^q^] 


where 


hr/i ± /3??3 ± 7V?1 
, 0772 ± fiys ± 7Vi 


u>2 = —aoR 


= a ± /? ± 7 — 


avi ± / 3 ? 7 | ± 777J 

{ar]2 ± /ills ± 


± Pvi ± 7ril 


(14) 


(15) 


In principle, the low-energy Hamiltonian reveals two 
possible scenarios, for the case of inhomogeneous arc¬ 
shaped strain which provides a Landau’s gauge for the 


H±=E± + ^[(a ± /3 ± 7){Qx + 9y) 

+ zf{x'^+ y'^)-2zf{xqy - yqx)] (16) 

.2 I 

where zf = ^[ar]2^l3'n3^7Vi] s-nd = mq( 0772 ±/3 773 ± 

7771 ). Notice that since z]*" < 0 in the case of M 0 S 2 , the 
energy spectrum obtained from (16) corresponds to the 
HO in the valence band and DQW in the conduction 
band, in a similar manner than in the case of the Landau 
gauge discussed before. It is, however, important to no¬ 
tice that the lack of the quadratic modulation, i.e. w^qf., 
of the HO band in Eq. (16) results in the appearance of 
the flat PLL in the valence band, as it has been dis¬ 
cussed by Cazalilla et al. [36]. This has to be compared 
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to the case of the arc-shaped deformation, for which such 
a quadratic term appears in (14) through and leads 
to a set of equidistant parabolic bands in the spectrum. 


A. Harmonic Oscillator solution in the valence 
band 


Using the single-band model for an arc-shaped defor¬ 
mation given by Eq. (14), it is possible to hnd an ap¬ 
proximated analytical expression for the energy disper¬ 
sion shown in Fig. 4. By using the basic physics of 
harmonic oscillator, it is straightforward to obtain the 
characteristic equidistant parabolic bands of this system, 
as given by 


Ev{n,qx) = E_- 



+ 


Amo 


Wo q 


2 

X 


^-7)1 



(17) 


as graphene case. This is so because the basis spinors 
in graphene are associated to the sublattice degree of 
freedom, while in monolayer TMDs they account for the 
conduction and valence band basis. The potential profile 
is shown in Fig. 5(a) for different values of qx- Notice 
that the potential is a symmetric (an asymmetric) double 
quantum well potential for qx = 0 {qx 7^ 0) with a barrier 
in the middle of the sample and two wells located at the 
edges. Therefore, the appearance of two wells close to 
the boundaries indicates the formation of a double quan¬ 
tum well (DQW) in the conduction band. According to 
Fig. 4, the energy dispersion becomes parabolic for en¬ 
ergies higher than a certain critical value at qx = 0. This 
feature obviously depends on the height of the parabolic 
barrier V{y,qx)- In fact, for energies higher than the 
barrier height, which is ~ 18.4 eV x {Ly/2R)'^, carrier 
motion does not be much affected by the existence of the 
barrier. 

At finite momentum qx the energy difference between 
two minima at y = AzLy/2 is given by 


where n = 0,1,2,.... Such a spectrum is consistent with 
the numerical result shown in Fig. 4. A simple estima¬ 
tion based on Eq. (17) suggests a separation between 
subbands of the order of AEy « 11.125^ eV. Since 
the maximum strain in the arc-shaped system is about 
£max ~ Ly/2R, therefore « 860 x ^ x . 

For a low temperature T = 5K and given system size 
Ly « 50 nm and strain strength about Smax = 0.1, we 
have AEy « 20 ^bT which implies sufficiently spaced 
levels as to be observed via STM spectroscopy. 


B. Double quantum well solution in the conduction 

band 

In order to understand the energy dispersion of the 
conduction band, we first include two hard walls at y = 
±Lyj2 which lead to the following Hamiltonian 

H+ = ^^{a + P + -/)ql + V{y,qx) (18) 

where 


SE = V{^,qx)-V{-^,qx) 


W2 Ly 

2 mo 


^ Qx 
( 20 ) 


which mimics an asymmetric DQW at any finite qx for 
which the two wells are no longer identical. 

If we neglect the hard-wall boundary condition, the 
Hamiltonian will be exactly solvable for the conduction 
band (18) which corresponds to the inverted harmonic 
oscillator equation [61, 62]. By performing elementary 
quantum mechanical approaches, we find the following 
eigenvalue relation for the wave function and correspond¬ 
ing energy eigenvalues in the conduction band 


CP In 
dz^ 4 


(j){z,e) 


e(j){z, e) 


( 21 ) 


where 

^ = Vw(y - w^qx) 

-1 \Amo{E - E+) + 2 

oj{a + P + 'y) h? ^ 


( 22 ) 


V{y,qx) =E+ + 

(19) 

where Q{x) is the step function, and Vq ^ 1 stands for 
the hard wall potential. Such hard-well boundary condi¬ 
tion can be realized by using external gates. Moreover, 
we expect that this boundary condition is justified for 
ribbons whose termination does not mix the valley de¬ 
gree of freedom, like the zigzag ribbon case. However, as 
it has been recently shown in Ref. 60, the boundary con¬ 
dition in the continuum model of TMDs is not as simple 


and uj = 2^—Wo/{oi +jd+ ^). This differential equa¬ 
tion is one of the standard differential equations for the 
parabolic cylindrical functions [63] whose two indepen¬ 
dent solutions can be written in terms of the confluent 
hypergeometric function, M (a, b, c) 

_ £i / .£ 11 

4’even(z, s) = 6 ^'^(~* 2 ”^ 4 ’ 2 ’*~^/ 

(f>odd(z,i) = ze-^^M -t 

which are even and odd functions with respect to the 
parity operator. The general solution of this prob¬ 
lem is a superposition of even and odd eigenfunctions. 
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(j) = ci(j)even + C 2 <^odd, where Ci ^2 are unknown constants. 
To find the corresponding eigenvalues, we must apply 
the hard-wall boundary condition which implies that the 
wave function must be zero at y = ±Lyj2. In this regard, 
it is easy to find the following relation for the eigenvalue 
problem 

ven (^2 ) 0 even ('\/^( ^2 ); ^) 

4>odd{V^{^ - w:^qx),e) <^odd(x/w(^ -I- w^qx),e) 

(24) 


At = 0 this eigenvalue problem reduces to search¬ 
ing for zeros of the confluent hypergeometric function. 
These solutions can be labeled as e = e„. If the 
even wave function satisfies the boundary condition 
4>even{'/^Ly/2,en) = 0, then we will get 


M 



1 1 

4’2’ 8 ) 


= 0 . 


(25) 


Otherwise, if the odd wave function becomes zero at the 
boundary, (poddiV^Ly/2, in) = 0, we will thus have 


M 



3 3 u}Ll\ 


= 0 . 


(26) 


Notice that the solutions of the above eigenvalue problem 
depend on both the microscopic details of the electronic 
band structure, which enter in the definition of the pa¬ 
rameter uj, and the ribbon width Ly in the form in{d}Ly). 
Once the set of solutions £„ are obtained, the correspond¬ 
ing energy levels at = 0 in the conduction band can 
be evaluated using Eq. (22). The energy levels are given 
by 

Ec{n,q^ =0)=E+- ^—yJ\w^{a + f3 + -f)\in{ujLl) 

(27) 

We numerically check that the lowest band has even sym¬ 
metry, as expected for a symmetric DQW potential. If we 
take qx 0 then the conduction band Hamiltonian will 
not commute with parity symmetry, since the potential 
is asymmetric for finite For any finite g^,, we solve Eq. 
(24) numerically for the lowest hve energy bands and re¬ 
sult is shown in Fig. 5(b) which is consistent with the full 
numerical calculation of the coupled differential equation 
(see Fig. 4). This agreement proves that our single-band 
picture, which predicts the HO and DQW physics in the 
valence and conduction bands, respectively, are appropri¬ 
ate models for the low-energy physics of strained TMDs. 
Furthermore, we emphasize that the theory presented 
here is general and can be straightforwardly adapted to 
other families of semiconducting 2D crystals. 


C. Effect of the scalar potential 

In this section, we consider the effect of the scalar po¬ 
tentials on the conduction and valence bands. We start 



y/ao 



FIG. 5. (Color online) (a) Potential profile in the condnction 
band, (b) Five low-energy bands around the K-point in the 
conduction band, as calculated from the single-band model 
(18) for the parameters Ly = 50ao and R = l.SLy. 


by performing full TB calculations in a zig-zag ribbon of 
monolayer M 0 S 2 and the results are shown in Fig. 6(a) 
and (b) for the unstrained and strained cases, respec¬ 
tively. These results are compared with those results 
obtained from the low-energy models given by Eq. (7). 
First, one notices the existence of three edge modes in the 
spectrum, in agreement with previous results [46, 57]. 
In Fig. 6 (b) we show the results for the band disper¬ 
sion in the arc-shaped strained case, which contains a set 
of roughly parabolic valence and conduction bands. In¬ 
terestingly, the energy gap around A'-point is no longer 
direct, which is one of the interesting consequences origi¬ 
nating from the scalar potential associated with this pro¬ 
file of the strain. Moreover, the crossing edge modes 
survive to the application of strain, whereas the flat high 
energy edge mode eventually enters into the bulk spec¬ 
trum for a higher strength of the strain. Intriguingly, 
one can see that the inter-level spacing in the conduction 
and valence bands for the strained system (Fig. 6(b)) in¬ 
creases dramatically as compared to the unstrained case 
(Fig. 6(a)), indicating that the origin of these levels does 
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not depend on the finite size effects but instead of the 
bulk potential induced by strain. In fact, they are the 
eigenvalues of the inverted and ordinary harmonic oscil¬ 
lator Hamiltonian of the electrons in the conduction and 
valence bands, respectively. Finally, Fig. 6(c) shows the 
results from the two-band low-energy model (7) for the 
same system, which are in good agreement with the cor¬ 
responding results from the full tight-binding model (Fig. 
6 (b)). 

In order to have some analytical insights of the effect 
of the scalar potential in the band structure of Fig. 6, 
we notice that such scalar potential for this strain profile 
can be written as 


D± {y) = pAtf y -k nfy'^] 

4mo 

4too 

Ogi? 

8 mo 

'"2 ^ li? op?^' 


(28) 


To include these potentials in the analytical calculation 
based on Eq. (14), we do need to replace wf with 
where 


V-^ — W-^ 

± 

^2 — ± 

V 3 =wi+WiW^'^-vfvf^. (29) 



Using the above relations, the energy dispersion in the 
valence band can be easily calculated. According to the 
negative sign of ap one can see that are negative 
for any value of strain, which means that the DQW and 
HO physics in the conduction and valence bands, as dis¬ 
cussed above, are still valid for TMDs with arc-shaped 
deformation in the presence of the scalar potential. 

Intriguingly, the shift of the conduction band minimum 
from the if-point should originate from the /qx term 
in V 2 according to the boundary condition equation given 
by Eq. (24). In particular, the band edge energy under 
arc-shaped strain can be expressed as 


^'VBM — 


Ag-A ^ ^ “o 


Amo 



/3-7) 



(30) 


The last term proportional to k)” is independent of strain 
strength. The strain-independent term should originate 
from our approximations. 


VI. STRAIN INDUCED VALLEY SHIFT IN 
HOMOGENEOUS DEFORMATIONS 

Using the strain induced modification of the hoppings 
given by Eq. (4), one can calculate the band dispersion 
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FIG. 6. (Color online) Energy dispersion in the presence of 
gauge fields and deformation potentials, (a) Full tight binding 
calculation in an unstrained zigzag ribbon with {Ly = 299ao). 
(b) Full tight binding calculation in a deformed zigzag ribbon 
with {Ly — 299ao , R ~ TLy). Dashed red line indicates a 
small shift of conduction band minimum with respect to the 
maximum of the valence band at the K-point. The energy 
gap around A-point is no longer direct originating from the 
scalar potential associated with this profile of the strain. The 
interlevel spacing in the conduction and valence bands for the 
strained system (panel b) dramatically increases as compared 
to the unstrained case (panel a), indicating that the origin 
of these levels does not depend on the finite-size effects, (c) 
Band structure calculated from the low-energy model (7) on a 
ribbon with hard wall boundary condition with {Ly = 299ao, 
R — 7Ly). Notice spin-orbit coupling has not been considered 
in this figure. 
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FIG. 7. (Color online). Energy dispersion of a uniaxially 
deformed monolayer M 0 S 2 calculated by using the full TB 
model. Solid (dashed) lines indicate spin up (down) compo¬ 
nents. The vertical dashed lines mark in addition the position 
of the conduction band minimum and valence band maximum 
according to the low-energy model (32)-(33). 

in the presence of different kinds of strain profiles. A 
teachful example is the case of a uniform uniaxial strain. 
It is commonly known, from DFT calculations [64], that 
uniaxial and shear strain induces a shift of the band edges 
from the K points, similar to the strained graphene. 
We can address this issue in our tight-binding approach. 
The energy dispersion of the monolayer M0S2, uniax¬ 
ially deformed along x, is shown in Fig. 7 close to 
the K- point [65]. Here, solid (dashed) bands indicate 
spin up (down) components. Both (conduction and va¬ 
lence) band edges shift in phase towards the F point for 
compressive {sq < 0 ) strain, whereas they move in the 
opposite direction for tensile strain, in agreement with 
DFT simulations [64]. A useful quantification of these 
valley-shifts with the strain can be obtained by means 
of an effective low-energy model. In the case of uniax¬ 
ial strain applied along x-direction, we have a^e/hAi = 
rjiAx ,aoe//iA 2 = r] 2 Ax, and aoe/hA^ = rf^Ax where 
A = {1 — v)£o, where v is the Poisson ratio. This uniax¬ 
ial strain shifts the valleys along the cc-direction, hence 
we safely set qy = 0. In this case, the approximated 
Hamiltonian for spin up around the AT-point is 

H = 4—^(1 ~ ^ 2 ) + I? + tokia^ 

+ ^ ^ {akl+Pa^kl), ( 31 ) 

dmoag 

where ki = q + rjiA and q = a^qx- Notice that just 
the leading term of the spin-orbit coupling is taken into 
account. The Hamiltonian can be easily diagonalized to 
obtain its band dispersion. Thus, it is straightforward 
to find the position of the conduction band minimum 


(^CBm) and the valence band maximum (^vem), which 
are given by 


ar]2 + I3r]3 -f 7771 

Qcbm — ^ , , 

a + p + ^ 

(32) 

arj2 - Pm - 7771 

^VBM ^ ' n ? 

a — p — 7 

(33) 


where the scalar potentials have no contribution to the 
leading term of the valley-shift. Importantly, in the par¬ 
ticle and hole bands, which they have the same effective 
mass (a = 0 ), gvBM = 9cbm = —A and the position 
of the valence and conduction band extreme are equally 
modified by strain. A similar behavior is obtained when 
m = m = Va = in which gvBM = gcBM = -qA. Since 
none of the previous special conditions apply to the case 
of strained M 0 S 2 , we expect different shifts for the elec¬ 
tron and hole band edges. Indeed, based on the numeri¬ 
cal value of the parameters in the low-energy model, we 
find gvBM = 0.76A and gcBM = 0.51A. Such a differ¬ 
ent strain induced a shift of the band edges leads to a 
direct-to-indirect gap transition in M 0 S 2 under uniaxial 
strain. 

It is interesting to compare the above results from the 
low-energy model with those results obtained from the 
TB. We do so in Fig. 7, where the vertical dashed lines in¬ 
dicate the position of the conduction band minimum and 
valence band maximum as obtained from the low-energy 
model (32)-(33). In the case of compressive strain, there 
is a good quantitative agreement between the two meth¬ 
ods. In the case of tensile strain, although the qualitative 
behavior is well captured by the low-energy model, the 
position of the valence band edge differs in the two cases. 
The reason is that, according to TB results, tensile strain 
enhances trigonal warping of the valence band, which is 
not considered within the simple low-energy model. A 
similar analysis can be done to understand valley-shifting 
induced by shear strain. In this case the result is similar 
to the one for uniaxial strain, with the difference that 
A = — 2eo and the deformation is along the y-direction. 
Finally, we remember that for biaxial strain, since A = 0, 
thus there is no valley-shift. 


VII. SPIN-STRAIN COUPLING 

Another interesting effect which is worth to be ad¬ 
dressed is the direct coupling between the spin and strain. 
Due to the spin-orbit coupling, we can manipulate spin 
degree of freedom of carriers just by controlling their or¬ 
bital motion. Using the direct spin-strain coupling, one 
can locally control the spin of carriers via the mechanical 
probe. In the absence of the spin-orbit coupling, there is 
no way to touch spin degree of freedom by deformation 
of the lattice via mechanical probe such as AFM tip. 

Our effective low-energy model (7) shows that the spin- 
orbit coupling terms are affected by external strain. At 
q = 0 the spin-orbit coupling in the conduction and va- 
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lence bands are modified as follows 

AA± = + Ag|774p ± A'|7?5p]|Alp + + toly) 

+ (34) 

These terms give rise a direct coupling between the spin 
and strain at the iiT-point. This coupling can allow for 
spin relaxation when translational symmetry in the sys¬ 
tem is broken because of, e.g., the existence of ripples 
which act as a long-range disorder potential. In fact, 
for the most studied case of graphene, such disorders 
lead to a spatially random spin-orbit coupling that might 
have implications for the spin relaxation [66]. The ef¬ 
fect of out-of-plane and in-plane deformation on the spin 
relaxation in the systems with similar symmetry have 
been studied [67, 68]. Although this effective spin-strain 
coupling seems to be weak, it is still comparable with 
some other energy scales like Zeeman energy or weak 
spin-splitting of the conduction band in semiconducting 
TMDs. This effect is expected to be especially relevant 
for the W families of TMDs (like WS 2 or WSe 2 ), for 
which the spin-orbit splitting of the valence band is more 
than twice the value for MoS 2 . From the experimental 
point of view, this tunable spin-orbit coupling via strain 
can be detected, in principle, using a photoluminescence 
measurement [69] because in the strained sample, we ex¬ 
pect a change in the relative position of A and B exciton 
peaks as compared with their position in the undeformed 
case. 

Finally, we notice that for finite value of q, another 
spin-dependent term 2 aoe//i{AQq.A 4 -|-A'q.A 5 (T 2 }®s ap¬ 
pears in the presence of strain. This term is important 
when A is finite and q is small enough (jA] > JaoqDj 
because in this regime we have [ttoqUA] > [aoq]^. There¬ 
fore, this new momentum dependent term in (7) originat¬ 
ing from strain (which is oc q) is dominant as compared 
to the (/-dependent term in the unstrained case, which is 
oc q^. 


VIII. SUMMARY 

In summary, we have studied the strain dependence of 
monolayer TMDs band structure, starting from a Slater- 
Koster tight-binding method which contains the neces¬ 
sary orbital contribution to describe the valence and 
conduction bands in the whole BZ. For a general in¬ 
homogeneous strain profile, we further calculate a low- 
energy Hamiltonian up to second order in momentum, 
strain and rotation tensors. Our numerical and ana¬ 
lytical calculations, based on TB and continuum mod¬ 
els, show a strong particle-hole asymmetry in the en¬ 
ergy spectrum of the system. We have shown that the 
momentum dependent terms in the low-energy model of 
monolayer TMDs acquire different strain induced vector 
potentials. According to our calculations using the low- 
energy model, the electronic spectrum of deformed sin¬ 
gle layer TMDs reveals a gauge-dependance which implies 


that these strain induced vector fields cannot be referred 
as gauge fields. Consequently, the simple pseudomag- 
netic field picture like, which is well-known in the case 
of strained graphene [37] is no longer valid in deformed 
monolayer TMDs. 

We have applied our theory to calculate the band struc¬ 
ture in the illustrative case of arc-shaped deformation. 
We show that this profile of strain induces three main 
fictitious gauge fields in Landau’s gauge. The dispersion 
relation of the system within this gauge shows several 
band crossings in the electron sector of the spectrum, 
and equidistant parabolic subbands in the hole spectrum. 
Our analytical calculations show that the energy spec¬ 
trum in the conduction and valence bands originate, re¬ 
spectively, from the solutions for a double quantum well 
and a harmonic oscillator Hamiltonian. This DQW and 
HO physics in two bands is also expected for a triangular 
strain profile. 

Finally, we study the shift of the conduction and va¬ 
lence band edges of MX 2 in the presence of homogeneous 
strain, finding a transition from direct to indirect gap. 
Moreover, the coupling between spin degrees of freedom 
and strain has been analyzed. We show that this ef¬ 
fect can be considered as a correction on the spin-orbit 
coupling that can useful for strain engineered spintronic 
applications. 
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Appendix A: On-site and hopping matrices 

In this Appendix, we provide the analytical expres¬ 
sions for the different contributions to the tight-binding 
Hamiltonian (3). The on-site terms of the Hamiltonian 
can be written in a compact form [10]: 

' = ("o" « ) ■ 
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where 


Ao 0 0 \ 

Em = j 0 62 —i^MSz 1 , 

yo iXmSz 62 J 

(^P+txx 0 \ 

ex = i^Sz ^p + tyy 0 . (A2) 

\ 0_0 e^ - tjj _ I 


Here, Am and Ax stand for the spin-orbit coupling of M 
(metal) and X (chalcogen) atoms, respectively [10] and 
= ± indicates z-component of spin degree of freedom. 
The terms = tyy = Ypp-jr, t'tz = ^ppa take into account 
the effects of the vertical hopping Vpp between the top 
and bottom chalcogen atoms. 

Below we list the hopping terms of the model. For the 
nearest neighbor hopping between M and X atoms we 
have [25] 


J.MX 

''1 


7^7 


/—9VpdTr + V^Vpda 

I 5y/3VpdTr + SVpda 

\—VpdTV — iV^Vpda 


SX^YpdTz Ypda 
9Vpd7r 

SVStpdTr + Y/pdrj 


Y2.VpdT^ + \/iVpdcr \ 

— 2'/3VpdTT + ^Ypda 1 

(^YpdTT ~ ^\/3Vpd(j / 


.MX 

^2 


.MX 

^3 


A. 

7V7 


7/7 


— 6^/3VpdTT + 2Vpdc + V^Vpda\ 

— QVpdTV — 4-\/3Hpd(T 4:V^VpdTV — QVpda 

0 0 / 


(= 

\lWpdTr 

/ 9VpdT^ — VSVpda- 

I —5'^VpdTr — iVpda 

V —VpdTT — ^V^Vpda 


Ypda 

9Ypd7r ^^y^Ypda 

7^^/3Ypd7r ^Ypdcr 


12Vpd7r + V^Vpda \ 

— 2\/WpdTr + iVpda 

— dVpdn + Sy/SVpdaJ 


(A3) 


(A4) 


(A5) 


Next nearest neighbor hoppings correspond to pro¬ 


cesses between the same kind of atoms, M-M or X-X 
(see Fig. 1), and they are given by 
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.MM 

^1 


1 

4 


^ SVddS + VddcT 

^{ — VddS + Vddcr) 

\ —^{VddS — Vdda) 


^{ — VddS + Vddcr) 

\{VddS + 12Vdd,r + SVddcr) 

^{VddS-V/dd.+iVdda) 


2 {VddS Vdda) \ 
^{VddS-^d.+Wdda) 
\{^VddS + 4:Vdd7r + ^Vdda) ) 


.MM 

^2 


1 

4 


/ “V/ddS + Vdda 
\/i{VddS — Vdda) 

V 0 


'/^{YddS — Vdda) 0 \ 

VddS + iVdda 0 

0 Wdd..) 


,MM 

^3 


1 

4 


^ "iVddS + Vdda 
^{ — VddS + Vdda) 
\ ^{VddS — Vdda) 


^{ — VddS + Vdda) 
\{VddS + ISVddTT + 3Vdda) 
-^{VddS - ^Vdd. + Wdda) 


2 {VddS Vdda) 
-^(VddS-Wddrr+SVdda) 
\(,^VddS + 4Vdd7r + ^Vdda) 


,XX 

^1 


1 

4 


( ^Vpprr “t“ Vppa 
'V^iVppT^ Vppu) 
0 


'n/S ( VppT^ Vppa ) 

Vpprr 4“ ‘^Vppa 
0 



(A 6 ) 


(A7) 


(A 8 ) 


(A9) 


^2 = 0 Vpp^ 0 

V 0 0 Vpp^j 


(AlO) 


,xx 

^3 


1 

4 


^Vpprz A Vppa 
\^{VppTY Vppa) 
0 


-\/3{VppTr - Vppa) 0 \ 

VppTY 4“ SVppfj 0 I 

0 Wpp^J 


(All) 


The direction of the hopping indicated by subindexes 1,2, 
and 3 can be seen in Fig. 1. 


Appendix B: Derivation of the low-energy 
Hamiltonian 

In this Appendix, we calculate the low-energy Hamilto¬ 
nian around the high symmetry K points using Lowding 
partitioning method [34]. The approach is similar to 
the one used in Ref. [42]. Here, the local strain is in¬ 
troduced by means of a local change of the two-center 
Slater-Koster matrix elements as a consequence of the 
local modulation of the interatomic bond lengths. We 
assume a general form of the inhomogeneous deformation 
with a large wavelength. To consider such deformation, 
we use the following relations for the bond lengths 
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r^^=ax 


= a\ 


rf^=ax 


1 Uxx ^yx \ / 


+ 


^yy ^xy 


2v^/ V2v^ 2^3 2 


+ 4 


_l_ ’^yy 


VsJ vVa V3 


1 ^xx ^yx \ I ^ ^ 

2 ”^ 2 ”^ 273 / ' 2^3 ' 2 


_l_ "^yy _|_ ^xy ^ _|_ 1 


MM(XX) _ 

V — Qj 


\ 


V3' 


Uy 


MM(XX) 


1 ^ 

2 2 2 

\2 


~Y 2 


yy ^xy 

Y 


1 + Uxx) + Ux 


xy 


MM(XX) _ 

Tg —a 


i 


1 Ux 


y/^U. 


yx 


I _|_ f _|_ _|_ Yy \ 


(B.l) 


in which Uij = duj/dxi. 

To obtain the low-energy model given in Eq. (7), we 
follow the next eight steps: 

1. We take the six-band tight-binding Hamiltonian (3) 
for a given spin subspace in the deformed system 
as follows 


iJ(k,u) = i7TB[6 X 6 ] (B.2) 


in which u is a tensor with matrix elements Uij. 

2. We expand the Hamiltonian (B.2) around the K 
point of the BZ, K = 47 r/ 3 a(l, 0 ), obtaining 


i7(k,u)«i7o + i7i(0 (B.3) 


where Hq = i7(K,0) and ^ = 

{qx,qy,Uxx,Uyy,Yy,Uyx}, where q = k - K 
and |q|a <C 1 . is obtained as 




dH(K + q, u) 


d^i 


5=0 


2 


.d^H{K + q, u) 


d^f 

d^H{K + q,u) 




C=o 


4=0 


(B.4) 


In this expansion, we assume that and com¬ 
mute, which is the case of the homogenous de¬ 
formation. The generalization for the inhomoge¬ 
neous strain case can be done, for long wavelength 
strain profiles, by replacing ujxy, and QiSjk with 
e,j{r), uJxy{r), and -i{dr;ejk{r) + ejk{r)dr-)/2, re¬ 
spectively. 


3. To find the low-energy Hamiltonian in the subspace 
of the conduction band minimum (CBM) and va¬ 
lence band maximum (VBM), we do a transforma¬ 
tion in orbital space using the unitary operator C/q, 
which diagonalizes Hq . After solving the eigenvalue 
problem 


Horn = Km (B.5) 

we obtain 

uo = im. m, m, m. m. m) ( b . 6 ) 

where we have ordered the eigenstates such that 
IV'o) and I^/^q) correspond to the lowest conduction 
and highest valence band for the given spin, respec¬ 
tively. Then, we apply this unitary transformation 
to move from orbital basis {H) to the band basis 
{H') as 


H'= U^o + HiiOPo. (B.7) 


4. We do the following replacement in H' 

^xx ^XXJ 

'^yy ~ 

^xy — ^xy T ^xyi 

Uyx — ^xy ^xy- (B. 8 ) 

Now, the Hamiltonian H' is a function of wave 
vector, q, symmetric strain tensor, e, and anti¬ 
symmetric rotation tensor, uj. 
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5. We decompose H' into two parts H' = Hd + V 


Hd = 


V = 


_fhn[4x4] 0 

0 /l22[2 X 2] 

0 /ri2[4 X 2] 

hl^[2x4] 0 


(B.9) 


where hn is defined in the subspace 
{IV'o), IV'o)’l'^o)> IV’o)}, whereas ^22 is defined 
in the subspace {|'0o)> IV'o)}- Notice that the high 
energy hn and low-energy /122 blocks in (B.9) is 
coupled with the off-diagonal element V. The 
analytical expression of the block components, i.e. 
hij, are too lengthy to be included here. 


We emphasize that the extension to the inhomoge¬ 
neous case is already done in (7) just by considering the 
local nature of the spin-independent contribution to the 
scalar potential (T>), and its spin-dependent part (<5A), 
as well as non-zero commutation of the momentum and 
the fictitious vector fields (i.e. [q, A^] ^ 0 ). This kind of 
extension from homogeneous to inhomogeneous deforma¬ 
tion is common in studying strained conventional semi¬ 
conductors and graphene [70, 71]. 

Note that there is also a trigonal warping term which 
is not included in Hamiltonian (7). For the unstrained 
case, it has the form 

Hu, = tmoQ • CT*tTa:q • cr* -k t 2 alT{ql - 3q,uqy){a' + 

(B.13) 


6 . An additional unitary rotation is performed to 
project V into each of these subspaces. We em¬ 
ploy the quasi-degenerate perturbation theory by 
using as rotation operator. This allows to 
drop the first-order V in the transformed Hamil¬ 
tonian, H" = e~^H'e^ = Hd + V + [Hd,0] + 
[V, O] + 5 [[Hd, 0],0] + - ■ ■, leading to the constraint 
V + [Hd, O] = 0. The generator of the transforma¬ 
tion takes the form 


O = 


f ° 

X 4] 


^[4 X 2 ]\ 

0 J 


(B.IO) 


where 77/122 — hiiV = /112 is solved to find the 77 
matrix as a function oi{qx,qy, e^x, Syy, exy,0Jxy} up 
to second order for the given spin index. 


where ti = —0.14eV, t 2 = leV, a’ = 0.44, (3' = —0.53. 
Here, the trigonal warping term contains three parame¬ 
ters (a',/3',/i). It is easy to show that all these terms 
combine with each other to lead the characteristic con¬ 
tribution z± cos 3(/) to the low-energy dispersion at the 
K-point, where z± = / 2 (a'±/3') ± 2/oti/(A — A_rs), and 
z+{z-) stands for the conduction (valence) band. 

It is interesting to notice that the direction of warping 
in both bands is opposite if z+Z- > 0 , and it is the 
same otherwise. If a' = 0 the warping in both bands are 
in the same direction and with same warping strength. 
Furthermore, a and a' are the sources of asymmetry in 
effective masses and trigonal warping directions between 
the conduction and the valence band, respectively. In 
our case, z+z- < 0 which means same warping direction 
in the two bands. 


7. Then, H" = Hd + ^[V,0] + • ■ ■ is our final ef¬ 
fective Hamiltonian with two decoupled subspaces. 
Following a straightforward calculation, the effec¬ 
tive Hamiltonian of the low-energy bands can be 
obtained for a given spin index as follows, 


H 2 b {qx 3 qy, ' 


xxj ^xyt '-^xy 


) = 


h 22 + ^{77^/112 + hl ^ r ]}. 


8 . Finally, we consider the relations 


(B.ll) 


IIe[Al] — ^xx ^yy, 
Im[Al] = -2exy, 


V — Sxx + e 


yy 


(B.12) 


and factorize (B.II) to reach 
of the Hamiltonian Eq. (7). 
simply extract the numerical value of 
parameters in our low-energy model i.e. 
{to, ti,t 2 ,a, /3, Ao, A, Ao, A, a', (3', Aq, A', 771 , 772 , 773 , 774 , 
af, a^, Qfg , 02 ^, from the numerical 

values of the pre-factors, and the result is given in 
Table H. 


the form 
Then, we 
all 


Appendix C: Derivation of the single band model 

Since Ai is weak as compared to the other vector po¬ 
tentials entering in our theory, we can apply perturba¬ 
tion methods to decouple the conduction and the va¬ 
lence bands, driving the two-bands Hamiltonian (7) into 
two one-band Hamiltonians as given in Eq. (14). To 
do this we use a canonical transformation similar to 
what we have done to get (7) in Appendix B, in which 
Hd = dia,g[hc,hv] and V 12 = 1 ^ 2 ^ = hcv By neglecting 
the spin splitting in the conduction band and the momen¬ 
tum dependence of the spin-splitting in the valence band, 
we obtain an expression equivalent to (B.9) by defining 
the following relations 


hr = 


An 


n? 


T-i I ^^0 2 I o 2 

D+ + -^^1^2 + 




V5, 


, Ao — A bar 2 bar „ 2 

hy — -h D- -I- A_ -|- — 'j^P'^3 

_ t^ ^ 

ibcv — ^ 1 


(C.l) 
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j-2 

where b = , K = hn^hy = / 122 , and hcv = ^ 12 - To 

calculate the generator of the transformation 


O = 



(C.2) 


we must first find d, which obeys the following relation 


'& =-hcv{hc-hy) ^ - ['d,hc]{hc-hy) ^ (C.3) 


''0 


A - A_ 


+ - ^(|q+ |A 3 |= - 7 |q+ |a.P) 

(C7) 

where 7 = 4mou^/(A — A_). We can also consider higher 
second order terms originating from the expansion (C.3), 
with the result 


ffp) = 


1 




(A-A_)2 ^~lg 

Pivi + V3){A^3 + TrqTTi)] - 2^7r|7r37ri} 


(C. 8 ) 


We solve Eq. (C.3) in an iterative perturbative approach 
in the low-momentum and low-strain limits. In this re¬ 
gard, it is easy to show that 


d = 


^( 1 ) ^( 2 ) 

A-A_ ^ (A-A_)2 


(C.4) 


where 


for the conduction band, and 

+ Pirn + V3){'^iT^l + ttsttI)] -H /3(7ri7r|7r3 -h 7r|7ri7r|)} 

(C.9) 

for the valence band, where we have neglected the con¬ 
tribution from the scalar potential term D±. 

Appendix D: Numerical eigenvalue problem 




(C.5) 


Therefore, performing the canonical transformation on 
the two-band Hamiltonian (7) we obtain, to first order in 
hcv, two decoupled Hamiltonians for the conduction and 
valence bands, Hy^y) = where 


To discritize and find the eigenvalues of the Hamil¬ 
tonian for an arc-shaped monolayer TMDs, we use the 
method of moments [72] with a trigonal basis to sat¬ 
isfy the hard wall boundary conditions which require 
vanishing of the wave function at the boundaries. We 
can expand the wave function in a set based on trigo¬ 
nal basis T„{y) = T{y - y„) as ipc = Y^n'^nTniy) and 
4>y = YnXnTniy), where (pc and (py refer to the conduc¬ 
tion and valence band spinor components respectively, 
and 




\y\ < Ly/iN + 1)-J^ 
otherwise ^ ' 


hW = 


Aq 


D. 


m 


A-x_il 


7 ^(a|q+ -h^|q-h -|-7|q-|- ^Aip) 

4777-0 n h fi 

(C. 6 ) 


and yn = Ly{jP^ ~ k)- ^o^e that Ly and N are the 
width of the system along the y-direction and the number 
of basis functions, respectively. It should be mentioned 
that using this trigonal basis guarantees the zero value 
of the wave function at the boundaries. To do this nu¬ 
merical calculation we need to know the following matrix 
elements in this set of trigonal basis 
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{Tm\dy\Tn) — —2{N + l)6m,n + (N + l){6m,n+l + ^m,n-l) 
{Tm\dy\Tn) = -^{Smn+1 ~ ^m,ra-l) 


(r™|T„) = 


2 1 
^m,n “t“ ^(^m,n+l “t“ 


(1 + iV) [3 ’ 6 


{Tm\y\Tn) — 


2n-l-N ^ ‘2.n-N ^ 2n-2- N ^ 

2 *^m,n—1 “1“ 1 


(7;n|2/2|T„) = 




(1 + 1V)3 


10(n-^)2 + l 20(n-f)2 + l 


15 


+ 


120 


^m,n+l “1“ 


20(n- + 1 

120 


^m,n—1 


(D.2) 


Using these matrix elements one can discretize the Hamil¬ 
tonian (11). Since trigonal basis has non-zero overlap, us¬ 


ing this method we obtain a generalized eigenvalue prob¬ 
lem that we solve numerically. The results can be seen 
in Figs. (3) and (4). 
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